今日やること!

7.1 アメリカ合衆国の州単位での地図

まずはプロットしたいデータの確認から

2016年のアメリカ大統領選挙の、州ごとの投票率および各得票数に関する変数のデータを可視化する。

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(socviz)

election %>% 
  select(fips, state, total_vote, r_points, pct_trump, party, census) %>% 
  slice_sample(n = 5)
## # A tibble: 5 × 7
##    fips state                total_vote r_points pct_trump party      census   
##   <dbl> <chr>                     <dbl>    <dbl>     <dbl> <chr>      <chr>    
## 1    17 Illinois                5594825  -16.9       38.4  Democratic Midwest  
## 2    11 District of Columbia     311268  -86.8        4.09 Democratic South    
## 3    42 Pennsylvania            6166729    0.710     48.2  Republican Northeast
## 4    26 Michigan                4824542    0.220     47.2  Republican Midwest  
## 5    48 Texas                   8993166    8.98      52.1  Republican South

FIPS (federal information processin standard) コードは、州・準州に割り当てられる連邦情報処理標準による番号。

party_colors <- c('#2E74C0', '#CB454A')
p0 <- ggplot(data = subset(election, st %nin% 'DC'),
             mapping = aes(x = r_points, 
                           y = reorder(state, r_points),
                           color = party)) +
  geom_vline(xintercept = 0, color = 'grey30') +
  geom_point(size = 2)
p0

p1 <- p0 +
  scale_color_manual(values = party_colors) +
  scale_x_continuous(breaks = c(-30, -20, -10, 0, 10, 20, 30, 40),
                     labels = c('30\n(Clinton)', '20', '10', '0',
                                '10', '20', '30', '40\n(Trump)')) +
  facet_wrap(~census, ncol = 1, scales = 'free_y') +
  guides(color = 'none') +
  labs(x = 'Point Margin', y = '') +
  theme(axis.text = element_text(size = 8)) +
  theme_bw()
p1

「空間データを必ずしも空間構造として表現しなくても良い」

地図上で描画しよう

library(maps)
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
us_states <- map_data('state')
class(us_states)
## [1] "data.frame"
str(us_states)
## 'data.frame':    15537 obs. of  6 variables:
##  $ long     : num  -87.5 -87.5 -87.5 -87.5 -87.6 ...
##  $ lat      : num  30.4 30.4 30.4 30.3 30.3 ...
##  $ group    : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ order    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ region   : chr  "alabama" "alabama" "alabama" "alabama" ...
##  $ subregion: chr  NA NA NA NA ...
head(us_states)
##        long      lat group order  region subregion
## 1 -87.46201 30.38968     1     1 alabama      <NA>
## 2 -87.48493 30.37249     1     2 alabama      <NA>
## 3 -87.52503 30.37249     1     3 alabama      <NA>
## 4 -87.53076 30.33239     1     4 alabama      <NA>
## 5 -87.57087 30.32665     1     5 alabama      <NA>
## 6 -87.58806 30.32665     1     6 alabama      <NA>
dim(us_states)
## [1] 15537     6

まずは白地図を書く。

p <- ggplot(data = us_states,
            mapping = aes(x = long, y = lat, 
                          group = group)) +
  geom_polygon(fill = 'white', color = 'black')
p

州ごとに色分けしたコロプレス図を描いてみる。

p <- ggplot(data = us_states,
            mapping = aes(x = long, y = lat, 
                          group = group,
                          fill = region)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  guides(fill = 'none')
p

デフォルトではメルカトル図法で書かれているが、広域のGISデータの描画では歪みが生じる。アルベルス正積円錐図法を使うことで、面積の比率が正しくなる。

p <- ggplot(data = us_states,
            mapping = aes(x = long, y = lat, 
                          group = group,
                          fill = region)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  coord_map(projection = 'albers', lat0 = 39, lat1 = 45) + # lat0, lat1はアメリカの場合のデフォルト
  guides(fill = 'none')
p

選挙のデータを、left_join()関数で地図データに結合する

先程の州ごとのデータに、州名をキーにして、選挙関連の統計情報をdplyr::left_join()関数で結合する。

head(election)
## # A tibble: 6 × 22
##   state      st     fips total_vote vote_margin winner party pct_margin r_points
##   <chr>      <chr> <dbl>      <dbl>       <dbl> <chr>  <chr>      <dbl>    <dbl>
## 1 Alabama    AL        1    2123372      588708 Trump  Repu…     0.277     27.7 
## 2 Alaska     AK        2     318608       46933 Trump  Repu…     0.147     14.7 
## 3 Arizona    AZ        4    2604657       91234 Trump  Repu…     0.035      3.5 
## 4 Arkansas   AR        5    1130635      304378 Trump  Repu…     0.269     26.9 
## 5 California CA        6   14237893     4269978 Clint… Demo…     0.300    -30.0 
## 6 Colorado   CO        8    2780247      136386 Clint… Demo…     0.0491    -4.91
## # … with 13 more variables: d_points <dbl>, pct_clinton <dbl>, pct_trump <dbl>,
## #   pct_johnson <dbl>, pct_other <dbl>, clinton_vote <dbl>, trump_vote <dbl>,
## #   johnson_vote <dbl>, other_vote <dbl>, ev_dem <dbl>, ev_rep <dbl>,
## #   ev_oth <dbl>, census <chr>
head(us_states)
##        long      lat group order  region subregion
## 1 -87.46201 30.38968     1     1 alabama      <NA>
## 2 -87.48493 30.37249     1     2 alabama      <NA>
## 3 -87.52503 30.37249     1     3 alabama      <NA>
## 4 -87.53076 30.33239     1     4 alabama      <NA>
## 5 -87.57087 30.32665     1     5 alabama      <NA>
## 6 -87.58806 30.32665     1     6 alabama      <NA>
# electionの州名を、us_statesと同様に全て小文字に変換する
# <-- よくある変換。そしてエラーが出やすい
election$region <- tolower(election$state)

# region列をキーにして、us_statesデータにelectionデータを結合する。
#   left_join(): 左のus_statesデータを優先して、electionデータを結合する
#   right_join(): 右のelectionデータを優先して、us_statesデータを結合する
#   inner_join(): us_statesとelectionデータに共通する列だけ結合する
#   full_join(): どちらか一方に含まれる列が結合される。存在しない列はNAになる
us_states_elec <- left_join(us_states, election, by = 'region')
head(us_states_elec)
##        long      lat group order  region subregion   state st fips total_vote
## 1 -87.46201 30.38968     1     1 alabama      <NA> Alabama AL    1    2123372
## 2 -87.48493 30.37249     1     2 alabama      <NA> Alabama AL    1    2123372
## 3 -87.52503 30.37249     1     3 alabama      <NA> Alabama AL    1    2123372
## 4 -87.53076 30.33239     1     4 alabama      <NA> Alabama AL    1    2123372
## 5 -87.57087 30.32665     1     5 alabama      <NA> Alabama AL    1    2123372
## 6 -87.58806 30.32665     1     6 alabama      <NA> Alabama AL    1    2123372
##   vote_margin winner      party pct_margin r_points d_points pct_clinton
## 1      588708  Trump Republican     0.2773    27.72   -27.72       34.36
## 2      588708  Trump Republican     0.2773    27.72   -27.72       34.36
## 3      588708  Trump Republican     0.2773    27.72   -27.72       34.36
## 4      588708  Trump Republican     0.2773    27.72   -27.72       34.36
## 5      588708  Trump Republican     0.2773    27.72   -27.72       34.36
## 6      588708  Trump Republican     0.2773    27.72   -27.72       34.36
##   pct_trump pct_johnson pct_other clinton_vote trump_vote johnson_vote
## 1     62.08        2.09      1.46       729547    1318255        44467
## 2     62.08        2.09      1.46       729547    1318255        44467
## 3     62.08        2.09      1.46       729547    1318255        44467
## 4     62.08        2.09      1.46       729547    1318255        44467
## 5     62.08        2.09      1.46       729547    1318255        44467
## 6     62.08        2.09      1.46       729547    1318255        44467
##   other_vote ev_dem ev_rep ev_oth census
## 1      31103      0      9      0  South
## 2      31103      0      9      0  South
## 3      31103      0      9      0  South
## 4      31103      0      9      0  South
## 5      31103      0      9      0  South
## 6      31103      0      9      0  South

結合済みのデータを描画しよう

選挙結果で共産党と民主党のどちらが勝利したのかを可視化してみる。

p <- ggplot(data = us_states_elec,
            aes(x = long, y = lat,
                group = group,
                fill = party)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  coord_map(projection = 'albers', lat0 = 39, lat1 = 45)
p

最初の作図の色と揃えてみる。ついでに見た目をcowplot::theme_map()で調整

library(cowplot)
p <- ggplot(data = us_states_elec,
            aes(x = long, y = lat,
                group = group,
                fill = party)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
  scale_fill_manual(values = party_colors) +
  labs(title = 'Election Results 2016',
       fill = NULL) +
  theme_map() +
  theme(legend.position = 'bottom')
p

ドナルド・トランプ氏の得票率を示してみる

p <- ggplot(data = us_states_elec,
            aes(x = long, y = lat,
                group = group,
                fill = pct_trump)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
  scale_fill_gradient(low = 'white',
                      high = '#CB454A') +
  labs(title = 'Trump vote',
       fill = 'Percent',
       fill = NULL) +
  theme_map() +
  theme(legend.position = 'bottom')
p

対立する2党への投票の結果を、中間点から分岐する配色で示してみる。

パターン1. 中間色を白にする。

p0 <- ggplot(data = us_states_elec,
            aes(x = long, y = lat,
                group = group,
                fill = d_points)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
  labs(title = 'Winning margins',
       fill = 'Percent') +
  theme_map() +
  theme(legend.position = 'bottom')

p1 <- p0 +
  scale_fill_gradient2()
p1

パターン2: 中間を紫

p2 <- p0 + 
  scale_fill_gradient2(low = party_colors[2],
                       mid = scales::muted('purple'),
                       high = party_colors[1],
                       breaks = c(-25, 0, 25, 50, 75))
p2

パターン3: 民主党が強い基盤を持っているワシントンD.C.を除外して可視化

p3 <- ggplot(data = filter(us_states_elec, region %nin% 'district of columbia'),
            aes(x = long, y = lat,
                group = group,
                fill = d_points)) +
  geom_polygon(color = 'grey90', size = 0.1) +
  coord_map(projection = 'albers', lat0 = 39, lat1 = 45) +
  scale_fill_gradient2(low = party_colors[2],
                       mid = scales::muted('purple'),
                       high = party_colors[1],
                       breaks = c(-25, 0, 25, 50, 75)) +
  labs(title = 'Winning margins',
       fill = 'Percent') +
  theme_map() +
  theme(legend.position = 'bottom')

p3

コロプレス図は、マッピングしている変数を部分的にしか表現できない。地図+地図に頼らない可視化はセットで作ると良い。

さて、自分のGISデータで可視化するにはどうしたらいい?

従来: 方法はほぼ選択の余地なし

  • ポリゴンデータ

    • spパッケージ
  • ラスタデータ

    • rasterパッケージ

現在: 特にラスタデータのパッケージが乱立中

  • ポリゴンデータ

    • sfパッケージ
  • ラスタデータ

    • terraパッケージ

    • starsパッケージ

芳賀の経験から・・・

「自分が分析で使いたいパッケージに合わせて、ポリゴン・ラスタデータを扱うパッケージを選ぶと良い。」

  • tidyverseパッケージと相性がいいのはsfパッケージ & starsパッケージ (同じ作者です)

  • terraパッケージは従来のrasterパッケージの作者が開発したもので、ほぼ同じ機能が実装されているし、関数も豊富。

  • ややこしい作業をしないならばsf + starsでよいかも。

  • もし高度なラスタ演算が必要ならば、terraパッケージでラスタ演算–>starsで描画、という流れが良いかもね。

「自分が分析で使いたいパッケージに合わせて」ってどういうこと?

  • 分析で使いたいパッケージの仕様書に書かれている関数のヘルプを読むのです。

    • 例えば、景観分析でよく使うlandscapemetricsパッケージのcalculate_lsm()関数は、argumentsとして、Raster Layer, Stack, Brick, SpatRaster (terra), stars, or a list of rasterLayersを取れるので、terraでもstarsでも、rasterでも大丈夫。

    • (https://cran.r-project.org/web/packages/landscapemetrics/index.html)

7.alpha 日本のGISデータを可視化してみようか

北摂の行政区域のデータをRで可視化してみる

library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
city_boundary <- sf::st_read(here::here('data', 
                                        'N03-20220101_27_GML', 
                                        'N03-22_27_220101.shp'))
## Reading layer `N03-22_27_220101' from data source 
##   `/Users/ch/Documents/GitHub/tm-data-viz/data/N03-20220101_27_GML/N03-22_27_220101.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 144 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 135.0913 ymin: 34.27182 xmax: 135.7466 ymax: 35.05129
## Geodetic CRS:  JGD2011
class(city_boundary)
## [1] "sf"         "data.frame"
summary(city_boundary)
##    N03_001            N03_002            N03_003            N03_004         
##  Length:144         Length:144         Length:144         Length:144        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##    N03_007                   geometry  
##  Length:144         MULTIPOLYGON :144  
##  Class :character   epsg:6668    :  0  
##  Mode  :character   +proj=long...:  0
glimpse(city_boundary)
## Rows: 144
## Columns: 6
## $ N03_001  <chr> "大阪府", "大阪府", "大阪府", "大阪府", "大阪府", "大阪府", "…
## $ N03_002  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ N03_003  <chr> NA, "大阪市", "大阪市", "大阪市", "大阪市", "大阪市", "大阪市…
## $ N03_004  <chr> NA, "都島区", "福島区", "此花区", "此花区", "此花区", "此花区…
## $ N03_007  <chr> NA, "27102", "27103", "27104", "27104", "27104", "27104", "27…
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((135.0918 34..., MULTIPOLYGON (((…

北摂の市町村のみ抽出してみる

hokusetsu <- c('豊中市', '池田市', '箕面市', '能勢町', '豊能町',
               '吹田市', '高槻市', '茨木市', '摂津市', '島本町')
hokusetsu_city_boundary <- city_boundary %>% 
  filter(N03_004 %in% hokusetsu)
hokusetsu_city_boundary
## Simple feature collection with 13 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 135.3302 ymin: 34.7307 xmax: 135.6827 ymax: 35.05129
## Geodetic CRS:  JGD2011
## First 10 features:
##    N03_001 N03_002 N03_003 N03_004 N03_007                       geometry
## 1   大阪府    <NA>    <NA>  豊中市   27203 MULTIPOLYGON (((135.4478 34...
## 2   大阪府    <NA>    <NA>  豊中市   27203 MULTIPOLYGON (((135.499 34....
## 3   大阪府    <NA>    <NA>  池田市   27204 MULTIPOLYGON (((135.4392 34...
## 4   大阪府    <NA>    <NA>  吹田市   27205 MULTIPOLYGON (((135.5068 34...
## 5   大阪府    <NA>    <NA>  高槻市   27207 MULTIPOLYGON (((135.5838 34...
## 6   大阪府    <NA>    <NA>  茨木市   27211 MULTIPOLYGON (((135.546 34....
## 7   大阪府    <NA>    <NA>  茨木市   27211 MULTIPOLYGON (((135.5301 34...
## 8   大阪府    <NA>    <NA>  箕面市   27220 MULTIPOLYGON (((135.4734 34...
## 9   大阪府    <NA>    <NA>  摂津市   27224 MULTIPOLYGON (((135.5555 34...
## 10  大阪府    <NA>    <NA>  摂津市   27224 MULTIPOLYGON (((135.5985 34...

ggplotで市町村名ごとに可視化してみる。日本語が含まれるので、https://github.com/Gedevan-Aleksizde/fontregisterer を参考にしてfontregistererパッケージをインストールしておく。

library(fontregisterer)
## You can see available font families by `names(quartzFonts())`
sans <- fontregisterer::get_standard_ja_fonts()['sans']
p_hokusetsu <- ggplot(data = hokusetsu_city_boundary,
                      mapping = aes(fill = N03_004)) +
  geom_sf() +
  scale_fill_grey() +
  theme_minimal() +
  labs(fill = '市町村') +
  theme(text = element_text(family = sans), 
        axis.text.x = element_text(angle = 45, hjust = 1))
p_hokusetsu

北摂周辺の土地利用ラスタデータをプロットしてみる。

まずはデータを読み込む。

library(stars)
## Loading required package: abind
lulc <- stars::read_stars(here::here('data', 'L03-b-14_5235', 'L03-b-14_5235.tif'))
lulc
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  L03-b-14_5235.tif 
##  50     :389732    
##  70     :107739    
##  10     : 57694    
##  110    : 30196    
##  100    : 18851    
##  160    :  8083    
##  (Other): 27705    
## dimension(s):
##   from  to  offset        delta  refsys point values x/y
## x    1 800     135      0.00125 JGD2000 FALSE   NULL [x]
## y    1 800 35.3333 -0.000833333 JGD2000 FALSE   NULL [y]

まずはプロットしてみる。

p_lulc <- ggplot() +
  geom_stars(data = lulc) +
  scale_fill_manual(values = c('yellow', 'beige', 'springgreen3', 'orange', 'red',
                               'grey30', 'grey70', 'brown', 'blue', 'lightyellow',
                               'lightblue', 'lightgreen', 'white'),
                    labels = c('田', 'その他の農用地', '森林', '荒地', '建物用地',
                               '道路', '鉄道', 'その他の用地', '河川地および湖沼',
                               '海浜', '海水域', 'ゴルフ場', '解析対象外')) +
  theme_minimal() +
  labs(fill = '土地利用') +
  theme(text = element_text(family = sans), 
        axis.text.x = element_text(angle = 45, hjust = 1))
p_lulc

北摂周辺の行政区域データと土地利用ラスタデータの重ね合わせ

まず、市町村界を、rasterデータと同じJGD2000に投影変換する

hokusetsu_city_boundary_proj <- sf::st_transform(hokusetsu_city_boundary, 
                                                 crs = sf::st_crs(lulc))
hokusetsu_city_boundary_proj
## Simple feature collection with 13 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 135.3302 ymin: 34.7307 xmax: 135.6827 ymax: 35.05129
## Geodetic CRS:  JGD2000
## First 10 features:
##    N03_001 N03_002 N03_003 N03_004 N03_007                       geometry
## 1   大阪府    <NA>    <NA>  豊中市   27203 MULTIPOLYGON (((135.4478 34...
## 2   大阪府    <NA>    <NA>  豊中市   27203 MULTIPOLYGON (((135.499 34....
## 3   大阪府    <NA>    <NA>  池田市   27204 MULTIPOLYGON (((135.4392 34...
## 4   大阪府    <NA>    <NA>  吹田市   27205 MULTIPOLYGON (((135.5068 34...
## 5   大阪府    <NA>    <NA>  高槻市   27207 MULTIPOLYGON (((135.5838 34...
## 6   大阪府    <NA>    <NA>  茨木市   27211 MULTIPOLYGON (((135.546 34....
## 7   大阪府    <NA>    <NA>  茨木市   27211 MULTIPOLYGON (((135.5301 34...
## 8   大阪府    <NA>    <NA>  箕面市   27220 MULTIPOLYGON (((135.4734 34...
## 9   大阪府    <NA>    <NA>  摂津市   27224 MULTIPOLYGON (((135.5555 34...
## 10  大阪府    <NA>    <NA>  摂津市   27224 MULTIPOLYGON (((135.5985 34...

北摂の範囲だけの土地利用データにするために、lulcをclipする

hokusetsu_lulc <- sf::st_crop(lulc, hokusetsu_city_boundary_proj)
# ちょっと可視化したいだけならこれで十分
plot(hokusetsu_lulc)

では、オーバーレイしてみよう。

p_hokusetsu_lulc <- ggplot() +
  # ラスタデータの描画
  geom_stars(data = hokusetsu_lulc) +
  scale_fill_manual(values = c('yellow', 'beige', 'springgreen3', 'orange', 'red',
                               'grey30', 'grey70', 'brown', 'blue', 'lightyellow',
                               'lightblue', 'lightgreen', 'white'),
                    labels = c('田', 'その他の農用地', '森林', '荒地', '建物用地',
                               '道路', '鉄道', 'その他の用地', '河川地および湖沼',
                               '海浜', '海水域', 'ゴルフ場', '解析対象外'),
                    na.value = NA) +
  theme_minimal() +
  # ベクタデータの描画
  geom_sf(data = hokusetsu_city_boundary,
          fill = NA, color = 'black', size = .5) +
  # プロットの調整
  theme_minimal() +
  labs(fill = '土地利用') +
  theme(text = element_text(family = sans), 
        axis.text.x = element_text(angle = 45, hjust = 1))
p_hokusetsu_lulc
## Warning: Removed 62938 rows containing missing values (geom_raster).

市町村ごとの土地利用の割合を計算してみる

ラスタデータの各グリッドに、hokusetsu_city_boundary_projの市町村名の情報を空間結合する。

lulc_joined <- st_join(hokusetsu_lulc, hokusetsu_city_boundary_proj)
## Warning in st_join.stars(hokusetsu_lulc, hokusetsu_city_boundary_proj): st_join
## found 1508 1-to-n matches, taking the first match for each of those

空間結合した市町村名で各グリッドをグループ分けして、市町村ごとの土地利用を集計する。

lulc_joined_df <- 
  # ラスタデータ変換し、
  # 1行に1グリッドの土地利用・市町村名が格納されたtibble形式のデータフレームを作成する
  as.data.frame(lulc_joined) %>% 
  tibble() %>% 
  # 土地利用が格納された列名をlulcに変更
  rename(lulc = L03.b.14_5235.tif) %>% 
  # 市町村名がNAではないグリッドだけ抽出
  filter(!is.na(N03_004) & !is.na(lulc)) %>% 
  # 市町村ごとに土地利用クラスを集計
  group_by(N03_004, lulc) %>% 
  summarise(N = n()) %>% 
  # 土地利用クラスのIDを名前に変換
  mutate(lulcname = case_when(lulc == 10 ~ '田', 
                              lulc == 20 ~ 'その他の農用地', 
                              lulc == 50 ~ '森林', 
                              lulc == 60 ~ '荒地', 
                              lulc == 70 ~ '建物用地',
                              lulc == 91 ~ '道路', 
                              lulc == 92 ~ '鉄道', 
                              lulc == 100 ~ 'その他の用地', 
                              lulc == 110 ~ '河川地および湖沼',
                              lulc == 140 ~ '海浜', 
                              lulc == 150 ~ '海水域', 
                              lulc == 160 ~ 'ゴルフ場', 
                              lulc == 255 ~ '解析対象外',
                              TRUE ~ as.character(NA)))
## `summarise()` has grouped output by 'N03_004'. You can override using the
## `.groups` argument.

さあ、プロットしようか。

p_lulc_bar <- ggplot(data = lulc_joined_df, 
                     mapping = aes(x = N, y = reorder(lulcname, N),
                                   fill = lulc)) +
  geom_bar(stat = 'identity') +
  scale_fill_manual(values = c('yellow', 'beige', 'springgreen3', 'orange', 'red',
                               'grey30', 'grey70', 'brown', 'blue', 'lightyellow',
                               'lightblue', 'lightgreen', 'white'),
                    labels = c('田', 'その他の農用地', '森林', '荒地', '建物用地',
                               '道路', '鉄道', 'その他の用地', '河川地および湖沼',
                               '海浜', '海水域', 'ゴルフ場', '解析対象外'),
                    na.value = NA) +
  theme_bw() +
  facet_wrap(~N03_004) +
  theme(text = element_text(family = sans)) +
  labs(x = 'グリッド数', y = '土地利用区分') +
  guides(fill = 'none')
p_lulc_bar

最後に・・・R関係で困った時には

  1. パッケージのヘルプを見に行く (“[package name] package”)
  1. チュートリアルとかムックを探す
  1. 素直にgoogle検索 (“[package name] package” AND “[error message]”)

  2. 日本国内なら TokyoR に参加してみる https://tokyor.connpass.com/

  3. @hagachi に聞いてみる

Enjoy your data analysis!

チュートリアル&あなた自身のデータでまずは手を動かしてみて、たくさんのエラーと最後の成功体験を楽しんでくださいな。